#include<../data_tools/usefull_routines.c>

particulas *particles;

int main(int argc, char *argv[])
{
  int i,ndisk,dummy1=0,dummy2=0;
  double *radio,rvir,r200;
  double Mtotal;
  size_t *ind;
  char *infile,outfiles[1000];
 
  infile = argv[1];
  r200 = atof(argv[2]);  
  rvir = atof(argv[3]);
  

  FILE *galaxyfile;
  galaxyfile = fopen(infile,"r");
  ndisk = contadorParticulas(infile);

  particles = (particulas *)malloc((size_t)ndisk*sizeof(particulas));
  if(particles == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  radio = (double *)malloc((size_t)ndisk*sizeof(double));
  if(radio == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  ind = (size_t *)malloc((size_t)ndisk*sizeof(size_t));
  if(ind == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }


  Mtotal = 0.0;
  for(i=0; i<ndisk; i++)
    {
      fscanf(galaxyfile,"%d %f %f %f %f %f %f %f %f\n",
	     &particles[i].id,
	     &particles[i].pos[0],&particles[i].pos[1],&particles[i].pos[2],
	     &particles[i].vel[0],&particles[i].vel[1],&particles[i].vel[2],
	     &particles[i].mass,
	     &particles[i].pot);
      radio[i] = sqrt(pow(particles[i].pos[0],2) + pow(particles[i].pos[1],2) + pow(particles[i].pos[2],2));
      Mtotal = Mtotal+particles[i].mass;
      ind[i] = i;
    }

  //Computing properties of the particles
  printf("Computing properties of the particles\n");
  particlesProperties(ndisk);  
  printf("Done.\n");
  printf("\n");

  /*
  //Printing galaxy with all properties
  printf("Printing galaxy with all properties in *.cmb_all.*\n");
  sprintf(namefiles,"%s.cmb",infile);
  printAscii(namefiles,1);
  printf("Done.\n");
  printf("\n");
  */
  gsl_sort_index(ind,radio,1,ndisk);
  
  printf("Done.\n");
  printf("\n");

  //Computing radial mass profile
  printf("Computing mass radial\n");
  FILE *masafile;
  sprintf(outfiles,"%s_radial_mass.dat",infile);
  masafile = fopen(outfiles,"w");
 
  for(i=0; i<ndisk; i++)
    {
      fprintf(masafile,"%e %e\n",radio[ind[i]],particles[ind[i]].mr);
    if((radio[ind[i]]>=rvir) && (dummy2==0))
	{
	  dummy2 = 1;
	  printf("\nM en Rvir(%e) : %e\n",rvir,particles[ind[i]].mr);
	}
      if((radio[ind[i]]>=r200) && (dummy1==0))
	{
	  dummy1 = 1;
	  printf("\nM en R200(%e) : %e\n",r200,particles[ind[i]].mr);
	}
    }
  
  printf("Done.\n");
  printf("\n");


  free(particles);
  free(radio);
  free(ind);
  fclose(galaxyfile);
  fclose(masafile);

  return 0;
}
